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Abstract. Extensive N-body simulations are among the key means for the study of numerous astrophysical 
and cosmological phenomena, so various schemes are developed for possibly higher accuracy computations. 

We demonstrate the principal possibility for revealing the evolution of a perturbed Hamiltonian system 
with an accuracy independent on time. The method is based on the Laplace transform and the derivation 
and analytical solution of an evolution equation in the phase space for the resolvent and using computer 
algebra. 

PACS. 98.52.-b Stellar systems 


1 Introduction 

The appearance of powerful computers has made the N-body simulations among efficient tools [1] for the study of 
various astrophysical and cosmological problems. The N-body gravitating systems cover a broad class of astrophysical 
problems, from the evolution of the Solar system as of a nearly integrable problem, up to the dynamics of star clusters, 
galaxies, and galaxy clusters as non-integrable problems. The cosmological simulations include from the evolution of 
primordial density perturbations up to the formation of dark matter galaxy halos, etc. 

The N-body simulation activities have been developed in two conventional directions: 

1. Numerical integration of N equations of motion in order to follow the evolution of the system. The main difficulty 
of these investigations based on standard iterative methods of numerical integration of differential equations (Runge- 
Kutta, etc.) is the inevitable storage of errors with the increasing of the number of steps. This leads to rapid loss of 
the reliability of the results derived, the faster, the larger is the number of equations (i.e. the number of particles) 
and/or the longer is the duration of the calculations. 

2. Investigation of the character of motion in the systems, which includes the finding out of integrals of motion, 
analysis of appearing stochastic regions, etc. The main feature of this direction is the application of methods of 
dynamical systems, e.g. of Kolmogorov-Arnold-Moser (KAM) theory, Lyapunov exponents, KS-entropy, etc P]. N- 
body gravitating systems, for example, do reveal variety of complex, chaotic properties which play a crucial role in 
their evolution and structure; see e.g. [31IU[51[S1[71[S] . 

Speaking of numerical investigations of dynamical systems, both, smooth and discrete, one may note that the 
second direction seems to be particularly fruitful. Beginning from the pioneer study by Fermi, Pasta and Ulam in the 
early 1950s, it has led to unexpected and profound results such as the strange attractors, the Feigenbaum sequence 
of bifurcations, etc. As compared to those remarkable results, the numerical (iterative) study of the evolution of 
Hamiltonian systems, due to the reasons mentioned, looks moderately successful and reliable. The latter is equally 
true for such an important class of Hamiltonian systems as the so-called nearly integrable systems P] 

H(/,i?,/3) = iLo(/)+/3ifi(/,^?), (1) 

where /, d are the action-angle coordinates and /3 is the small parameter defining the non-integrable part of the 
Hamiltonian Hi. 

The well-known Henon-Heiles system P] modeling the stellar motion in a spiral galaxy is of this type; the same 
class of Hamiltonians can be attributed also to various aspects of planetary dynamics. The unique role of the system 
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(1.1) in the Hamiltonian mechanics is determined undoubtedly by the proof of the KAM theorem, according to which, 
when distinct conditions are satisfied, there exist invariant tori, 

/ = /o + w(t?,/3), 

t? = 'do +a;t+ u(d,/3), 

where u,v are periodic functions [5]. 

It is well known that the KAM theorem with its roots goes back to the main problem of the celestial mechanics, the 
stability and evolution of the Solar system. Although the KAM theorem does not directly reveal the fate of the Solar 
system, its ideas were used for efficient numerical studies of the Solar system dynamics; see jTlIH] and other studies by 
Laskar et al. The KAM theorem itself does not determine the value of the perturbation /3, for which its statements are 
true, and regarding the planetary dynamics there are also principal difficulties in checking of the necessary conditions 
of the theorem with respect to each consideration. Moreover, due to the Arnold diffusion - a universal instability 
peculiar to non-linear systems of dimension N > 2 - irrespective of whether the Solar system satisfies the conditions of 
the KAM theorem or not, it still cannot remain stable. Finding itself after an accidental perturbation in the stochastic 
region of phase space, the system can remain in it for an infinitely long time and therefore, the observed picture has 
to be destroyed anyway - the planets must either fall onto the Sun or fly away. So, those fundamental results at least 
do predict that, strictly speaking, the Solar system cannot last forever anyway, although the time scale of the latter 
instability (Arnold diffusion), according to estimates, far exceeds the Hubble time. 

The example described above shows, on the one hand, the universality of systems with a perturbed Hamiltonian; 
on the other hand, the difficulties in direct application of the KAM theorem for revealing of the evolution of real 
physical and astrophysical systems. 

Connecting certain hopes with computers and speaking of their potentialities, one must mention the importance 
of computer algebra methods, which, in our opinion, offer new prospects for the investigation of complex dynamics. 

Below, we will show that computer algebraic methods along with those of dynamical systems enable the principal 
possibility of investigation of the evolution of a system without the effect of the storage of errors, i.e. of an accuracy 
independent on time. The search of error-free numerical integration schemes can be considered among the key problems 
of stellar dynamics. Problem 5 in [ID]; see also mi. 


( 2 ) 

( 3 ) 


2 Evolution in the phase space 


The convenience of the investigation of a Hamiltonian system with small perturbation for action-angle variables is 
connected with the following. As is well known, the Hamiltonian of an n-dimensional integrable system with a phase 
space = {p, g} , e.g. in the case of free oscillators with perturbation 


H{p, q) 


n 


E 



(wW' 

2 


+ PU{q) , 


( 4 ) 


i.e. a system having n first integrals of motion in involution, enables one to transit from the variables (p, q) to that of 
action-angle variables (/, ’d) 

qk = (24/a;'=)i/2cosi?'=, 

. (5) 

Pfc = (24w'=)i/2sini?4 


Then, for /3 = 0 the equations of motion 


4 


dH 


w'=(/) + 


dH 
dik ’ 


( 6 ) 


have a trivial solution. 


4 = const , 

-I-= const. 


( 7 ) 


As it follows from the Liouville theorem, in this case the phase space of the system splits into compact manifolds 
diffeomorphic to n-tori 

Tor" = {i?* mod 27r}. (8) 

When /3 ^ 0, the Hamiltonian equations read 


fk _ dH{I,-d,P) 

^ , 

hk _ dH{I,«,l3) 

^ ~ dIk 


( 9 ) 
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where 


H{I,^,P) = Y,uj^h + PV{I,d). 


k=l 


However, instead of solving or analyzing these equations as is done conventionally, we will proceed as follows. 
One can observe that if (see [2], Section 39) 




then for any smooth function f[x) the following holds: 


Therefore, 


or 


where 


’ dx^ dx'^ ^ ' 


+ LaUx) = 0 

La = -^A\x)^^. 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


Considering the function /(/,i?) in the phase space {!,'&) and applying the above observation we find that the 
equation describing the evolution of this function at (/, i?) = x changing according to (6), has the form 


where 


*^/(a;, t) + L{j3)f{x, t) = 0, 


L{P)=Lo+PB, 

T _ _ A dH d 

~ ^ dik ’ 

B{V) = -i 


dV _d _JL 

dIk dIk 


By means of the Laplace transform 


/(a;. A) = / dtf{x,t)e 


iXt 


(16) 


(17) 


(18) 


one can rewrite equation (16) as follows: 


- A f{x, X)+L f{x, A) = -i f{x, 0) 


(19) 


/(a:, A) = iRxi/3)g{x ), 
g{x) = f{x,0 ), 


where the resolvent R\{13) is equal to 


RxiP) = {x-m) 


-1 


One can show that the following expression is fulhlled 


Rx{P)g{x) = (X-Lq- f3B) ^ g{x) 

= {1 - PRxB)-^ Rxg{x) 

= Eto akm[RxB]'^Rxg{x) + o{/3E , 

i.e. at fixed /3, by variation of we can reach a given accuracy. 

If 

g{I,d)= gk{I)e^^^A>^ 

kGZ’^ 

where 

< kjd >= ki'd’' , 


( 20 ) 

( 21 ) 


( 22 ) 


(23) 
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then 

Rx9iI,^) = {X-Lor^ 9(1,^) 

= {X + i < uj,d0 >)~^ 9(1, t?) 

(24) 

= Efcez-(A + i < a;, a, >)-' e*<M> 

= 9k{I) x-<k,ui>' 

i.e. (18) can he calculated relatively easily. 

We choose the function 9 {x) in the form 

- n n 

9(1, ^) = - 2 + E • (25) 

fc=i fc=i 

This choice is conditioned by two factors: first, it has a single maximum on i?" x Tor” and hence, its evolution has 
the same property; second, it allows an elementary representation in the form of a Fourier expansion. 

Having R\ 9 , by means of the inverse Laplace transform we get the desired function, 

fix, *)= e"*^‘/(a;, X) = - Rx{^) 9 ix ), (26) 

where the contour of integration is 


C = {a + is, —oo < a < oo, s = a > 0} . (27) 

Therefore, the evolution of the initial function of f{x,0) = 9 {x), is found. Estimating after this the maximum of 
the function f{x,t) = for each t we find the time evolution of /, d, thus determining the evolution of the 

Hamiltonian system for those time instances. 

The principal difference between this and the standard iteration methods of integration of Hamiltonian systems is 
clear; here the time t is a parameter of f{x,t) and does not influence the accuracy of calculation of the latter. 


3 Computations via the new scheme 

Firstly, let us note the importance of application of the computer algebra in our calculations. Due to the absence 
of errors at an integration and differentiation operation, the final results practically do not contain any errors. To 
demonstrate this, we consider a dynamical system given by a two-dimensional Hamiltonian of two oscillators with 
small perturbations, 

H = ojili-\-UJ 2 I 2 + piicos'di. (28) 

The initial function /(/, 'd, 0) was chosen in the form of (25). Then the evolution of that function by means of the 
procedure described was obtained, 

/(/,!?, t) = exp(ii?i)(/i - /i(0))/i -b iexp(-idi(0))[exp(-iwit) - 1] 

-|-exp(ii9i)(/i — Ji(0))/i -b \ exp(3zz?i — — idi(0))x 

x[exp(iwit) - A(/i -/i( 0 ))^ - A(/2 - 12(0))^ ^ ^ 

-b| exp(idi — iojit — idi(O)) -b ^ exp(id2 — iuj2t — *^2(0)). 

The evolution of the surface determined by the function at different instances of time and for Ii = 

const, I 2 = const is shown in Fig. 1. In Figs. 2 and 3 the variations of and Ii by the time for different initial 
conditions and values of j3 are shown. The trajectories are regular with 27r period. During computations, at the 
increase (decrease) of /3 by an order of magnitude, the calculation time varied proportional, at the accuracy 10“^. 
The performed computations were little time consuming (for i7, 2600 3.4 GHz processor of 6 GB memory); we plan 
to apply the approach to real physical systems and provide more results in forthcoming papers. 

Let us stress again that the accuracy of the computations of the trajectory does not depend on time, since no 
iterations are involved. 
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t = 2 
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Fig. 1. The time evolution of the surface determined by Eq.(29) for 7i = const, I 2 = const. 


4 Conclusions 

We demonstrated a principal way of finding out of the evolution of a dynamical system independent on time accu¬ 
racy; namely, the phase space point /i(0),/2(0), t?i(0), t?2(0) describing the initial state of the Hamiltonian system is 
transferred to an arbitrary t, without any storage of errors. The key feature of this method is that we transfer the 
function /(/, d, 0) and, therefore, solve an equation which is sufficiently easier (linear), as compared to the Hamiltonian 
equations. Moreover, the former equation is solved analytically and not by an iteration procedure. 

Another remarkable advantage of this procedure is the following: the initial state of a non-linear system during 
evolution can find itself in the stochastic region of phase space. The description of this is impossible by integration 
(even numerical) of the Hamiltonian equations, while our method allows us to investigate even this phenomenon. 
This aspect is particularly important for the study of N-body gravitating systems in view of their well-known chaotic 
properties. Then this method can open new principal possibilities for cosmological and extensive astrophysical N-body 
simulations. 
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Fig. 2. The time variations of i9i for initial conditions /i(0) = 1.0,/2(0) = 2.0 and i?i(0) = 3.0,i92(0) = 2.0 and several values 

oip. 







